home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / MATH1 / FITPOL.PAS < prev    next >
Pascal/Delphi Source File  |  1985-04-03  |  4KB  |  148 lines
  1. program fitpol;        { -> 295 }
  2. { Pascal program to perform a linear least-squares fit }
  3. { to the ratio of 2 polynomials }
  4. { with Gauss-Jordan routine }
  5. { Sperate modules needed:
  6.             GAUSSJ}
  7.  
  8.  
  9. const    maxr    = 20;        { data prints }
  10.     maxc    = 4;        { polynomial terms }
  11.  
  12. type
  13.     ary    = array[1..maxr] of real;
  14.     arys    = array[1..maxc] of real;
  15.     ary2    = array[1..maxr,1..maxc] of real;
  16.     ary2s    = array[1..maxc,1..maxc] of real;
  17.  
  18. var
  19.     x,y,y_calc    : ary;
  20.     resid        : ary;
  21.     coef,sig    : arys;
  22.     nrow,ncol    : integer;
  23.     correl_coef    : real;
  24.  
  25. external procedure cls;
  26.  
  27. procedure get_data(var x: ary;        { independant variable }
  28.            var y: ary;        { dependant variable }
  29.            var nrow: integer);    { length of vectors }
  30. { get values for n and arrays x,y }
  31.  
  32. var    i    : integer;
  33.  
  34. begin
  35.   { clausing factors }
  36.   nrow:=10;
  37.   x[1]:=0.1;    y[1]:=0.9524;
  38.   x[2]:=0.2;    y[2]:=0.9092;
  39.   x[3]:=0.5;    y[3]:=0.8013;
  40.   x[4]:=1.0;    y[4]:=0.6720;
  41.   x[5]:=1.2;    y[5]:=0.6322;
  42.   x[6]:=1.5;    y[6]:=0.5815;
  43.   x[7]:=2.0;    y[7]:=0.5142;
  44.   x[8]:=3.0;    y[8]:=0.4201;
  45.   x[9]:=4.0;    y[9]:=0.3566;
  46.   x[10]:=6.0;    y[10]:=0.2755;
  47. end;        { procedure get data }
  48.  
  49. procedure write_data;
  50. { print out the answers }
  51. var    i    : integer;
  52. begin
  53.   writeln;
  54.   writeln;
  55.   writeln('  I      X      Y      YCALC      RESID');
  56.   for i:=1 to nrow do
  57.     writeln(i:3,x[i]:8:1,y[i]:9:4,y_calc[i]:9:4,resid[i]:9:4);
  58.   writeln; writeln(' Coefficients errors ');
  59.   writeln(coef[1]:8:5,' ',sig[1],' constant term');
  60.   for i:=2 to ncol do
  61.     writeln(coef[i]:8:5,' ',sig[i]);        { other terms }
  62.   writeln;
  63.   writeln('Correlation coefficient is ',correl_coef:8:5)
  64. end;        { write_data }
  65.  
  66. {procedure square(x: ary2;
  67.          y: ary;
  68.          var a: ary2s;
  69.          var g: arys;
  70.      nrow,ncol: integer);}
  71. { matrix multiplication routine }
  72. { a= transpose x times x }
  73. { g= y times x }
  74. {$I C:SQUARE.LIB }
  75.  
  76. {external procedure gaussj(var b:    ary2s;
  77.                   y:    arys;
  78.               var coef:    arys;
  79.               ncol:        integer;
  80.               var error:    boolean);
  81. }
  82. {$I GAUSSJ.LIB }
  83.  
  84. procedure linfit(x,        { independant variable }
  85.          y: ary;    { dependent variable }
  86.          var y_calc: ary;    { calculated dep. variable }
  87.          var resid:  ary;    { array of residuals }
  88.          var coef:   arys;    { coefficients }
  89.          var sig:    arys;    { error on coefficients }
  90.          nrow:       integer;    { length of array }
  91.          var ncol:   integer);    { number of terms }
  92.  
  93. { least squares fit to nrow sets of x and y pairs of points }
  94. { Seperate procedures needed:
  95.     SQUARE -> form square coefficient matrix
  96.     GAUSSJ -> Gauss-Jordan elimination }
  97.  
  98. var    xmatr        : ary2;        { data matrix }
  99.     a        : ary2s;    { coefficient matrix }
  100.     g        : arys;        { constant vector }
  101.     error        : boolean;
  102.     i,j,nm        : integer;
  103.     xi,yi,yc,srs,see,
  104.     sum_y,sum_y2    : real;
  105.  
  106. begin        { procedure linfit }
  107.   ncol:=4;    { number of terms }
  108.   for i:=1 to nrow do
  109.     begin        { setup matrix }
  110.       xi:=x[i];
  111.       yi:=y[i];
  112.       xmatr[i,1]:=1.0;    { first column }
  113.       xmatr[i,2]:=-xi*yi;    { second column }
  114.       xmatr[i,3]:=xi;    { third column }
  115.       xmatr[i,4]:=-sqr(xi)*yi
  116.     end;
  117.   square(xmatr,y,a,g,nrow,ncol);
  118.   gaussj(a,g,coef,ncol,error);
  119.   sum_y:=0.0;
  120.   sum_y2:=0.0;
  121.   srs:=0.0;
  122.   for i:=1 to nrow do
  123.     begin
  124.       xi:=x[i];
  125.       yi:=y[i];
  126.       yc:=coef[1]+(-coef[2]*yi+coef[3]-coef[4]*xi*yi)*xi;
  127.       y_calc[i]:=yc;
  128.       resid[i]:=yc-yi;
  129.       srs:=srs+sqr(resid[i]);
  130.       sum_y:=sum_y+yi;
  131.       sum_y2:=sum_y2+yi*yi
  132.     end;
  133.   correl_coef:=sqrt(1.0-srs/(sum_y2-sqr(sum_y)/nrow));
  134.   if nrow=ncol then nm:=1
  135.   else nm:=nrow-ncol;
  136.   see:=sqrt(srs/nm);
  137.   for i:=1 to ncol do        { errors on solution }
  138.     sig[i]:=see*sqrt(a[i,i])
  139. end;    { linfit }
  140.  
  141.  
  142. begin        { main program }
  143.   cls;
  144.   get_data(x,y,nrow);
  145.   linfit(x,y,y_calc,resid,coef,sig,nrow,ncol);
  146.   write_data
  147. end.
  148.